clear
invP = invParams;

%Mesh
nc = 8;
invP.mesh = meshDM(10*ones(nc,1),10*ones(nc,1),10*ones(nc,1));
invP.activeCells = invP.mesh.CC(:,3)<sum(invP.mesh.hz)/2;

%Tx/Rx locations
invP.loc = [sum(invP.mesh.hx)/2 sum(invP.mesh.hy)/2 sum(invP.mesh.hz)/2+5];

%Time stepping
invP.disT = [50 1e-2];
   
%Primary Fields   
invP.h0 = 10000*calcH0_VMD(invP);

%Interpolation matrix
invP.Q = getInterpmat(invP);

%Sigma and tau models
invP.sigmaModel = 1e-8*ones(invP.mesh.nc,1);
invP.sigmaModel(invP.activeCells) = 1e-3;

invP.tauModel = zeros(invP.mesh.nc,1);
invP.tauModel(invP.activeCells) = 0.01;

%True model
etaTrue = zeros(invP.mesh.nc,1);
etaTrue(invP.activeCells) = 0.1;

dTrue = fwd_TDIP(invP,etaTrue);
d0 = fwd_TDIP(invP,zeros(invP.mesh.nc,1));

figure(1);
loglog(abs(dTrue),'-o');
hold on
plot(abs(d0),'-rx');
hold off

WmtWm = invP.mesh.DIV*invP.mesh.DIV'+0.00000001*speye(invP.mesh.nc);
uncert = 0.1*abs(dTrue) + 1e-10;

%%

m = zeros(invP.mesh.nc,1);
m_ref = zeros(invP.mesh.nc,1);
beta = 1;

[hs_pred,~,j_pred] = calcHs_TDIP(invP,m,false);
dpred = projHs(invP,hs_pred);

gp1 = (dTrue - dpred)./(uncert.^2);
gp1 = -JTvec_TDIP(invP,m,j_pred,gp1);
gp2 = beta*WmtWm*(m-m_ref);

g = gp1+gp2;

phi_d_old = norm((dTrue-dpred)./uncert);
phi_m_old = (m-m_ref)'*WmtWm*(m-m_ref);
phi_old = phi_d_old + beta*phi_m_old;

%%

alpha = 8.5e-4;
dm = -g;
m_guess = m+alpha*dm;

d_guess = fwd_TDIP(invP,m_guess);

phi_d_guess = norm((dTrue-d_guess)./uncert);
phi_m_guess = (m_guess-m_ref)'*WmtWm*(m_guess-m_ref);
phi_guess = phi_d_guess + beta*phi_m_guess;



figure(1);
loglog(abs(dTrue),'-o');
hold on
plot(abs(d_guess),'-rx');
hold off

%%
rhs = -g;
rhs(~invP.activeCells) = 0;

F = @(x)Hvec(invP,m,j_pred,60000*beta,WmtWm,x);
[dm,flag,relres,iter,resvec] = pcg(F,-g,1e-80,10,WmtWm);




m_guess = m+dm;

d_guess = fwd_TDIP(invP,m_guess);

phi_d_guess = norm((dTrue-d_guess)./uncert);
phi_m_guess = (m_guess-m_ref)'*WmtWm*(m_guess-m_ref);
phi_guess = phi_d_guess + beta*phi_m_guess;



figure(1);
loglog(abs(dTrue),'-o');
hold on
plot(abs(d_guess),'-rx');
hold off

figure(2);
plot(m_guess);


